#ifndef PARAMETER_H
#define PARAMETER_H
#include "./Grid.h"
#include "./Constants.h"
#include "./UnitConvert.h"
#include "./input.h"
#include "./Util.h"
#include "./Absord.h"

using std::ofstream;
using std::string;

namespace Parameter{
     //! Time-------------------------------------------
    const double t0 = UnitConvert::TimeAU(input::Time_0);
    const double T_pulse1_SI = input::wavelength1 / Constants::cc;
    const double dt_SI  = T_pulse1_SI / (input::SP_pulse1 * 1.0);
    const double dt = UnitConvert::TimeAU(dt_SI);
    const int    nt = input::SP_pulse1 * input::NC_pulse1;
    const double t_max = dt * nt;
    const int    fft_num = (nt%2 == 0)? int(nt/2)+1 : int((nt-1)/2)+1;

    //! Radial-----------------------------------------
    const double r0     = UnitConvert::LengthAU(input::Radius_0);
    const double dr_min = UnitConvert::LengthAU(input::Delta_r_min);
    const double dr_max = UnitConvert::LengthAU(input::Delta_r_max);
    
    //! Longitudinal-----------------------------------
    const double dz_SI = input::z_factor * dt_SI * Constants::cc;
    const double z0 = UnitConvert::LengthAU(input::z_0);
    const double dz = UnitConvert::LengthAU(dz_SI);
    const int    nz = round(input::distance / dz_SI);

    //! input laser1------------------------------------
    const double phase0 = input::phase0; //cep
    const double phase1 = input::phase1; //cep
    const double phase2 = input::phase2; //cep
    const double fwhm    = UnitConvert::TimeAU(input::FWHM);

	const int    index_omega1 = input::NC_pulse1; // 傅里叶变换后，入射场中心频率omega1在频率array中的索引
    const double Waist1  = UnitConvert::LengthAU(input::waist1);
    const double Focus1  = UnitConvert::LengthAU(input::focus1);
    const double E_peak1 = UnitConvert::ElectricAU(input::laser_intensity);
    const double Wavelength1 = UnitConvert::LengthAU(input::wavelength1);
    const int    nr1     = 1.5 * Util::Get_Nradius(z0,Waist1,Focus1,Wavelength1,dr_min,dr_max);

    //! input laser2--------------------------------------
    const int    index_omega2 = input::NC_pulse2; // 傅里叶变换后，入射场中心频率omega2在频率array中的索引
    const double Waist2  = UnitConvert::LengthAU(input::waist2);
    const double Focus2  = UnitConvert::LengthAU(input::focus2);
    const double E_peak2 = input::E_factor * E_peak1;

    const double Wavelength2 = UnitConvert::LengthAU(input::wavelength2);
    const int    nr2 = (E_peak2 == 0)? 0 : 1.5 * Util::Get_Nradius(z0,Waist2,Focus2,Wavelength2,dr_min,dr_max);
    const int    nr  = (nr1 >= nr2)? nr1 : nr2;
    int nr_out = nr;

    //! nonlinear parameter-----------------------------------
    double Kerr = UnitConvert::KerrAU(input::eta2);
    double Ng = UnitConvert::InitialDensityAU(input::density_of_neutrals);

    //! results.out-------------------------------------------
    ofstream filename_grid(input::output_folder + "grid.txt");
    std::string filename_Erz_accuracy = input::output_folder + "Erz_accuracy";
    
    ofstream intensity   (input::output_folder + "intensity.txt");
    ofstream intensity_r0(input::output_folder + "intensity_r0.txt");
    
    ofstream THzIntensity   (input::output_folder + "THzIntensity.txt");
    ofstream THzIntensity_r0(input::output_folder + "THzIntensity_r0.txt");
    
    ofstream plasmaDensity   (input::output_folder + "plasmaDensity.txt");
    ofstream plasmaDensity_r0(input::output_folder + "plasmaDensity_r0.txt");
    
    ofstream out_phase1(input::output_folder    + "phase1.txt");
    ofstream out_phase1_r0(input::output_folder + "phase1_r0.txt");
    ofstream out_phase2(input::output_folder    + "phase2.txt");
    ofstream out_phase2_r0(input::output_folder + "phase2y_r0.txt");
    
    ofstream Etz(input::output_folder + "Etz.txt");
    ofstream THzSpectrumyOut(input::output_folder + "THzSpectrum.txt");
    
    ofstream Erza(input::output_folder + "Erza.txt" );
	ofstream Erza1(input::output_folder + "Erza1.txt" );
    ofstream Erza2(input::output_folder + "Erza2.txt");

    ofstream Erza_r0(input::output_folder + "Erza_r0.txt");
    ofstream Erza1_r0(input::output_folder + "Erza1_r0.txt");
    ofstream Erza2_r0(input::output_folder + "Erza2_r0");

    //! fileout_parameter-------------------------------------
    void fileout_parameter(){
        std::ofstream parameter(input::output_folder + "parameter.txt");
        if(parameter){
            parameter << "laser_intensity=" << input::laser_intensity <<endl;
            if(E_peak2 ==0){
                parameter << "[Single Pulse:]"<< endl;
                parameter << "wavelength1=" <<std::fixed <<std::setprecision(0)
                          << input::wavelength1<<"nm"<<endl;
		        parameter << "index_omega1=" << index_omega1 <<endl;
                parameter << "focus1=" <<std::fixed <<std::setprecision(0)
                          << input::focus1 << "nm" << endl;

            }
            if(E_peak2!=0){
                parameter << "[Dual Pulse:]"<< endl;
                parameter << "wavelength1=" <<std::fixed <<std::setprecision(0)
                          << input::wavelength1<<"nm"<<endl;
		        parameter << "index_omega1=" << index_omega1 <<endl;
                parameter << "focus1=" <<std::fixed <<std::setprecision(0)
                          << input::focus1 << "nm" << endl;
                
                parameter << "wavelength2=" <<std::fixed <<std::setprecision(0)
                          << input::wavelength2<<"nm"<<endl;
		        parameter << "index_omega2=" << index_omega2 <<endl;
                parameter << "focus2=" <<std::fixed <<std::setprecision(0)
                          << input::focus2 << "nm" << endl;
		    }
            parameter << "distance=" <<std::fixed <<std::setprecision(0)
                      << input::distance << "nm" << endl;
                      
            parameter <<"----------------------------"<< endl;
            parameter << "nr_out=" << nr_out << endl;
            parameter << "Nraius=" << nr << endl;
            parameter << "Ntime=" << nt << endl;
            parameter << "Nz=" << nz << endl;
            
            parameter << "dt=" << std::fixed << std::setprecision(4)
                      << dt << "a.u."<< " " << dt_SI << " fs" << endl;
            parameter << "dr_min=" << dr_min << "a.u."
                      << " " <<std::fixed <<std::setprecision(0)
                      << UnitConvert::LengthSI(dr_min) << " nm" << endl;
            parameter << "dr_max=" << dr_max << "a.u."
                      << " " <<std::fixed <<std::setprecision(0)
                      << UnitConvert::LengthSI(dr_max) << " nm" << endl;
            parameter << "dz=" << dz << "a.u."
                      << " " <<std::fixed <<std::setprecision(0)
                      << UnitConvert::LengthSI(dz) << " nm" << endl;

            parameter << "T_max=" << t_max << "a.u."
                      << " " << std::fixed << std::setprecision(4)
                      << UnitConvert::TimeSI(t_max) << " fs" << endl;
        }
    }
}

#endif